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f-^ ■ Abstract. The quasicontinuum approximation [2j| is a method to reduce the atomistic degrees of 

freedom of a crystalline solid by piecewise linear interpolation from representative atoms that are 
nodes for a finite element triangulation. In regions of the crystal with a highly nonuniform deforma- 
ns) ' tion such as around defects, every atom must be a representative atom to obtain sufficient accuracy, 
but the mesh can be coarsened away from such regions to remove atomistic degrees of freedom while 
retaining sufficient accuracy. We present an error estimator and a related adaptive mesh refinement 
algorithm for the quasicontinuum approximation of a generalized Frenkel-Kontorova model that 
enables a quantity of interest to be efficiently computed to a predetermined accuracy. 
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1. Introduction 



The solution of the equations for mechanical equilibria of a crystalline solid modeled by a classical 
atomistic potential requires the computation of the interaction of each atom with all of the other 
I> ' atoms in its sphere of influence. Due to the high computational complexity, it is generally not 

possible to obtain numerical solutions for systems that are large enough to simulate long-range 
elastic effects, even for short-ranged potentials. However, the local environment of nearby atoms 
is almost identical up to translation, except in the neighborhood of defects such as cracks and 
dislocations. The quasicontinuum method utilizes this slow variation of the strain away from 
defects to approximate the full systems of equations of mechanical equilibrium by equations of 
^ ■ equilibrium at a reduced set of representative atoms [T], |3, lla, \2^ . 

More precisely, the positions of the full set of atoms are obtained by piecewise linear interpolation 
from representative atoms that are nodes for a finite element triangulation. A quasicontinuum 
energy is defined as a function of the positions of the representative atoms. The quasicontinuum 
Jh I method is then used to obtain a solution to a desired accuracy with a significant reduction in the 

computational degrees of freedom by coarsening the finite element mesh away from the singularities. 
Near defects, sufficient accuracy can only be obtained if all atoms are representative atoms. In 
contrast to continuum models, the mesh cannot be refined past the atomistic scale. 

Even higher efficiency is achieved by approximating the total energy of the atoms in a coarsened 
triangle by the product of the strain energy density and the area of the triangle (or its higher 
dimensional analogue). The strain energy density is obtained from the energy per atom in a lattice 
which is a uniform strain of the infinite reference lattice. The uniform strain in turn is determined 
from the displacement of the nodes of the respective triangle. The quasicontinuum approximation 
we obtain this way allows the coupling between a region that is computed as in fully atomistic 
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simulations and a region that is computed using the methods of piecewise Hnear finite element 
continuum mechanics. 

It is well-known that an efficient adaptive algorithm is highly dependent on the quantity of 
interest or goal of the computation, see [l|, l5(. In this paper, we adopt this goal-oriented approach 
to obtain the approximation of a quantity of interest to within a desired tolerance. 

The reliable and efficient utilization of the quasicontinuum method requires both a strategy to 
determine the decomposition of the fully atomistic system into atomistic and continuum regions 
and a strategy for refinement within the continuum region. In [3) y] , we have developed a goal- 
oriented error estimator and a corresponding adaptive algorithm to decide between the atomistic 
model and the continuum model. In this paper, we develop an error estimator and a corresponding 
adaptive algorithm for mesh refinement in the continuum region. 

We have chosen initially to investigate as a model problem a generalization of the classical 
Frenkel-Kontorova model. The potential energy includes next-nearest-neighbor interactions in ad- 
dition to nearest-neighbor interactions so that the continuum energy of a representative atom is 
different from the atomistic energy of a representative atom for a nonuniform strain. 

We note that in contrast to error estimators and adaptive algorithms for mesh refinement in 
classical continuum mechanics, our continuum region is coupled to the atomistic region. This is 
a considerably more complex setting than in purely continuum models with classical boundary 
conditions. Also, our mesh refinement algorithm restricts the representative atoms which serve as 
the mesh points to the sites of the atoms in the reference lattice, although this could be relaxed 
away from the atomistic-continuum interfacial region. 

Algorithms for adaptive mesh refinement for the quasicontinuum method have been proposed and 



investigated in numerical experiments for several mechanics problems in [Ill.ll7l.l23l|. An a posteriori 



error indicator for a global norm was analyzed and tested for a variant of our one-dimensional 



quasicontinuum method in 20(]. The goal-oriented approach to adaptive mesh refinement for the 
quasicontinuum method was first investigated in [ISL |l9| . Mathematical analyses of several variants 
of the quasicontinuum method have been given in [o, M, llO|, ll3, llj, l2]| . We refer to [J, lD, llJ, 122, 1251] 



for alternative atomistic-continuum coupling methods. 

2. Quasicontinuum Approximation 

In this section, we introduce the atomistic model and its quasicontinuum approximation. We 
refer to [4] for a more detailed description. 

We consider a one-dimensional system of 2M atoms whose positions are denoted by y = 
{y^M+i: ■ ■ ■ iVm) G M^^^. The potential energy of the atomistic system is described by a func- 
tion 

We split the energy into atom-wise contributions by means of 

M 

^"(y)= E ^"(y)- (2.2) 

i=-M+l 



In this paper, we consider a Prenkel-Kontorova type model [15|] which serves as a one-dimensional 
description of a dislocation. We expect that the a posteriori error estimators that we introduce and 
the corresponding adaptive refinement algorithms will be applicable to more general quasicontinuum 
models. For the Frenkel-Kontorova model, the atom- wise contribution, £f, consists of two parts, 

ff = 81^^ + £:f'"^. (2.3) 



ADAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 



• • • • 



-6a^ -Sag -4ag -3a„ -2^^ -a„ 



• • 



^0 2ag Sa^ 4ag Sa^ Sa^ 



misfit 
energy 

modeled 
layer 



substrate 



Figure 1: Frenkel-Kontorova model. The wells depict the misfit energy, S^ 



The elastic part, S^'^ , describes nearest-neighbor (NN) interactions and next-nearest-neighbor 
(NNN) interactions, whereas £^ ' models the misfit energy of a slip plane sitting on an unde- 
formed substrate, see Figure [TJ The two parts are defined as 



+ \k2{yi - yi-2 - 2ao)^ + \k2{yi+2 -Vi- 2ao)' 



^t (y) 



\kQ {yi - {i- l)aof , i 



-M + 1,...,0, 
1,...,M. 



(2.4) 



Here oq € M denotes the equilibrium distance, and the moduli ko, ki, and k2 describe the strength of 
the misfit energy and the elastic interactions, respectively. To ensure coercivity, we require /cq > 
and ki + 2/c2 > 2|A;2|. 

Next, the continuum energy 



for each atom i is derived from the atomistic energy, which leads to 

^i'^iy) = \h2{yi - Vi-i - aof + \ki2{yi+i - yi - a^f, 
^\ko{yi-{i-l)aof, i = -M + 1, . . . ,0, 



(y) 



\kQ {yi 



[i - Ijao) 

\2 



(2.5) 



(2.6) 



laQ) 



1,...,M, 



where ki2 = ki + 4/c2, see J2!] for the details. We note that £^'^{y) = £^'^{y) if yi+2 — Vi+i = 
Vi+i -yi = yi- Vi-i = Vi-i - yi-2, but E^^'^iy) ^ ^f'^(y) more generally. 

For each atom i, we decide whether this atom is modeled atomistically or as continuum. An a 
posteriori error estimator for this task has been derived in [3| • Let 



1 if atom i is modeled atomistically, 
if atom i is modeled as continuum. 



and 



l-5f. 



We define the atomistic- continuum energy to be 

M M 

f-(y):= Y. ^'i^t{y)+ E ^'i^i^) 



(2.7) 



(2.8) 



i=-Af+l 



-Af+1 
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The next step is to coarsen out unnecessary atoms within the continuum region. This way, we 
restrict the system to the remaining atoms, cahed the representative atoms, or briefly repatoms. 
An a posteriori error estimator to determine the optimal coarsening will be developed in this paper. 
Let ij for j = — A^+ 1, . . . , A^ be the index of the j-th repatom, where 2N is the number of repatoms. 
We require that 

l-N+l < ^-N+2 < ■ ■ ■ < ij < ij+1 < ■ ■ ■ < iN-1 < ^N (2.9) 

and 

i-N+i = -M + l, e^N+2 = -M + 2, eN-i = M-l, £n = M. (2.10) 

Then 

Uj := ij+i - ij (2.11) 

gives the number of atomistic intervals between the repatoms j and j + 1- We denote the vector of 
repatoms by y^'^. 

The quasicontinuum energy is obtained by implicitly reconstructing the missing atoms from the 
repatoms by piecewise linear interpolation, and then computing the atomistic-continuum energy 
from the reconstructed vector. The piecewise linear interpolation can be written as the matrix 
multiplication /"Sy^c^ where 

7-ag ^3 ~ ^ jaq _^ (O ^r>\ 

for j = — A^ + 1, . . . ,N and k = 0, . . . jVj, and /"^ = otherwise. Hence, the quasicontinuum energy 
is given by 

£i^{yi^):=g^^(^piyi^). (2.13) 

Summation formulas for an efficient computation of £'i^[y'i^) without having to explicitly recon- 
struct the non-repatoms have been derived in [2]. 

We describe boundary conditions at the two leftmost atoms and the two rightmost atoms. We 
take into account two instead of one atom at each end of the chain because of the NNN interactions. 
For the fully atomistic system and the atomistic-continuum system, the boundary conditions read 
as 

„,a _ ..be a _ be a _ be a _ be 

y-M+i — vii-: y-M+2 — yi2^ yM-i — yr2^ yM — Vri^ (':>^^\ 

ae be ae be ac be ae be \ ' J 

y-M+1 = yii , y-M+2 = yi2 > y&i-i = yr2, yn = yri, 

for given values yf{^, yj*!, y^2> ^^'^ Vri^ ^'^d similarly the boundary conditions for the quasicontinuum 
system are given by 

„,9C _ „,be qe _ be ae _ he qc _ he /n -, r.\ 

y-N+1 — yiii y-N+2 — yi2? yw-i — yr2^ yN — yn- [z.io) 

Next, we introduce the solution spaces 

y":=M2M, Vo^:=R^^'^-\ V^:=R^^, V^' ■.= R^^-^ (2.16) 

and the boundary operators 

J" : Vq^ -^ V and J" : V^ -^ V^ (2.17) 

by 



J^. ■= dij for i = -M + 1, . . . , M, j = -M + 3, . . . , M - 2, 
Jlj:=5ij for i = -N + l,...,N, j = -N + 3, . . . ,N - 2. 



(2.18) 
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Figure 2: Solution spaces and their operators for the fuUy-atomistic level (a) and the quasicontinuum level 
(q). 



V^ and Vq are the suitable spaces for the atomistic system and the atomistic-continuum system, 
both with and without boundary values, whereas V^ and Vq are the suitable spaces for the qua- 
sicontinuum system. J" and J"^ extend vectors from Vq and Vq, respectively, by zero boundary 
values. Note that the interpolation operator /"'^ maps V to V"". The spaces and their operators 
are illustrated in Figure [2j 

To implement the boundary conditions, we define the vectors 

bca ._ jaq bcq ^ ya 

We seek minimizers of the three different potential energies subject to the boundary conditions 
described above, that is, vectors 



(2.19) 



,bca 



y« e rVo^ + y'^'^'^ C y", y'^" G rV^^ + y"™ C y^ and f^^ G J'^V^ + y'''^^ C V^ (2.20) 



bcq 



which minimize the potential energies £"", S""^, and £'^^, respectively. If we decompose 



y'^ 



JV" + y' 



bca 



-ac 



J V + y' 



bca 



and 



yqc^jq^qc^y 



bcq 



(2.21) 



then the fully atomistic solution y"^ G Vq, the atomistic-continuum solution y'^'^ G Vq, and the 
quasicontinuum solution y'^'^ G Vq are characterized by 



rQC 



argmin£:''(JV + y^""), 

yeVo" 

argminf'^^(J"y + y*'"') 

yGK)" 

argminf«"(J'?y + y*'=«) 



(2.22) 



Next, we write the energies in matrix notation as 

^''(y) = |(y - Si'^fD^E'^Diy - a") + 1 (y - h-fK{y - b'^), 
S'^'iy) = Uy- Si'^^D'^E'^'Diy - a'^) + i(y - h^)'^K(y - b") 



[y) = i(y- Si'^fD'^E'^'^Diy - a'^) + 1 (y - b'^)^K(y 



(2.23) 

f9'=(y) = l{piy - ai''fD'^E''^D{P'^y - a") + i(rV - h''fK{P'iy - b'^). 

We refer to [4] for the precise and lengthy definitions of the respective matrices. In matrix notation, 
the vectors y"^, y"'^, and y'^'^ are the solutions of the linear systems 

M^y" = r, 
^ac^ac ^ pc^ ^2.24) 

^gcy9c ^ ^qc^ 
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where the matrices M"', M'^^, and M'^^ are given by 

M" := J'^^iD'^E'^D + K)J'', 

^ac ._ ja^{^JJ,T^acj^ ^ ^ya^ (2.25) 

M"^ := J'1^I'"i^{D'^E''^D + K)I'"iJ'i, 
and where the right-hand sides f", i""^, and i''^ are defined as 

-J'^^D'^E'^'^Diy^'"' - a") - J"^i^(y''='^ - b"^), (2.26) 

_jg^ J<^q^ ^Tj^acj-^/jaq bcq _ ^a\ _ jq^ jaq^ j^ ^ jaq bcq _ j^a-v 



cac 



3. Goal-Oriented Error Estimation 

T 

We estimate the error y^ — J^ /^-^ j^y^c -^^ terms of a user-definable goal function 

Q:V^^ M, (3.1) 

that is, we aim at estimating 

IQ(y"'=)-Q(^"^/"W)|. (3.2) 

We assume that Q is hnear. Hence there exists some vector q G Vq such that 

Q(y) = q^y Vy e V^. (3.3) 

We decompose the error (|3.3p as follows: 

IQ(y") - Q(J"^/"VV')I = IQ(y" - y'^') + Q(y"' - J'^^/^V^y^^)! 

< IQ(y" - y"') I + IQ(y"' - J"^rv«y'"=) | (3.4) 

where 



^a-ac ._ a _ ac 
^ac-qc ._ ac _ J°-^ pi Jly'i'^ 



(3.5) 



The first term |Q(e"~"'^)| constitutes the modeling error and has been treated in [2]. The second 
term \Q{e"''^~'^'^)\ describes the error due to mesh coarsening and will be treated in the following. 
To facilitate the error analysis, we define the dual problems 

T T (3.6) 



We then have the basic dual identity for the goal-oriented error 

= g"^^^ M"^ (y'^'^ - J"^ I"' J'^y"'^) (3.7) 

= g'^^^7^"'=(J'^^rVy«") 
with the primal residual given by 

7e«^(y) := pc - M^^y. (3.8) 
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However, this quantity is too expensive to compute. We cannot solve for the dual solution g"*^ on 
the atomistic scale, since this has the same computational complexity as solving for the uncoarsened 
solution y'*'^. To overcome this obstacle, we replace g'^'^ by a dual solution from a coarser space. 

T 

A first idea would be to use J" J'*'? j^g^c instead of g"'^. However, this turns out to be useless 
due to Galerkin orthogonality: 

Lemma 3.1 (Garlerkin Orthogonality). We have 
or equivalently 

T T 

Proof. Multiplying the equation for y""^ from (I2.24|) by J'^ /"'^ J"" from the left gives 

jqTjaqT ja^ac^ac ^ jq^ jaq^ ja^ac _ ^^^^y^ 

It is easy to see that 

jajoTjaqjq ^ jaq jq (3_^2) 

Applying this to the equation for y''^ from (j2.24|) leads to 

jqT jaqT ja^jacja^ jaqjq^qc ^ pc ^ ^^ ^^^ 

Subtracting (IXTHI) from (f3lT]) . substituting the definitions (f2:26|) of i"" and f^ and using (fXT2]l 
and p.igp gives 



jq-r-^aq-r ja^ac^^ac_ja 


Tjaqjq^qc^ 


_ jq"^ jaq'^ ja^ac _ ^qc 


= J^"^/"*^ J^^J"^ 


_DT^acjj^ybc _ ^a^ _ ^^^bc 


T T 

- J" r« 

= 0, 
which completes the proof. 


-D'^E^^DiP'^y''^'^ - a") - K{I 










(3.14) 



D 

Hence, replacing g"'^ in (13. 7p by J*^ jo-Q ji gif^ always gives a zero estimate for the goal-oriented 
error. We need to use the dual solution from some space which is finer than Vq to get a non-zero 
estimate for the goal-oriented error, but which is coarser than Vq to make it computable. 

To this end, we introduce an additional level of refinement and denote it as the partial continuum 
(pc) level. Similarly to the repatoms on the qc-level, we chose a set of 2N pc-level repatoms which 
is a subset of all atoms and a superset of the qc-repatoms. This means we choose indices ij for 
J = -N + 1,...,N such that 

t^_JV+l ^ ^—N+2 < ■ ■ ■ < ^J" < "]+l < • • • < -tjv-l ^ ^N (3.15) 

and 

i.N+i = -M + 1, £_^+2 = -M + 2, ^^„i = M - 1, lfj = M. (3.16) 

To ensure that the pc-level is actually a refinement of the qc-level, we require every qc-level repatom 
to be a pc-level repatom. Hence there exist indices /ij such that 

^i=4,, i = -A^ + l,...,iV. (3.17) 
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Figure 3: Levels of refinement: atomistic-continuum (ac), partial continuum (pc), and quasicontinuum (qc) 
with indices (above the respective chain) and intervals (below the respective chain). 



(3.19) 



Similar to the definition of Vj, we denote the number of atomistic intervals between two pc-level 
repatoms by 

%:=fj+i-lj, j=-N + l,...,N-l. (3.18) 

See Figure [3] for an illustration of the three levels of refinement and the corresponding variables. 
We define the pc-level solution spaces 

VP := M2^ and V^ := R^^-\ 

the interpolation operators 



> foi j = -N + 1,...,N -1, k = 0, 
otherwise, 



rap 

ij+kj 


i^j— k 


jap 

ij+k,3+l 


k 


jap 


= 


jpq 

Hj+k,j 


_ vj-l^^+k+l^,, 


^3 


jpq 




tPI 


= 



(3.20) 



> ioT j = -N + l,...,N -1, k = 0,...,jij+i- fij, 
otherwise, 



the restriction operator 



^Z-=^HJ for j = -N + l,...,N, j=-N + l 



, . . . 5 X T , 



N, 



and the boundary operator 

Jlj-=^iJ for i = -N + l,...,N, j = -N + 3,...,N-2. 
Note that the interpolation operators are defined in such a way that 1°'^ factors as 

jaq _ japjpq^ 



(3.21) 
(3.22) 
(3.23) 
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Figure 4-' Solution spaces and their operators for the ac-level, pc-level, and qc-level. 



Altogether, all solution spaces and the corresponding operators are depicted in Figure [H 
We solve the dual equation on the newly introduced partial continuum level: 

jVfPCgpc ^ jp'^pp^^q^^ (3.25) 

where 

^pc ._ jp'^pp'^^D^E'^^D + K)PPJP. (3.26) 

We use gP'^ to get an approximation of the goal-oriented error: 

Qiy""") - Q{r'^ p" J'^y"^) ~ gp"'^ jp'^ pp'^ m'"'{r^ p''jiy'"'). (3.27) 

Due to Galerkin Orthogonality (Lemma 13. ip . we can subtract any vector in V^^ from the right 
hand side expression without changing its value. Thus, we subtract from gP^ the linear interpolation 
of the nodal values of g^'^ at the qc-repatoms, JP IP'^R'^PJPgP^, to get a vector in Vq that is zero 
at the qc-repatoms. This leads us to the error estimator 



T] :-- 



[g'"' 
[g"' 



jp' piRQpjPgP'^] ' JP' pi' pTv^ir njiy"") 
jp'^ipiRipjpgp^f[jp'^pi'^rr^ - Mpyp"] 



(3.28) 



as an approximation to ((37 

In Section [5l we will construct a numerical method based on this error estimator. We will then 
compute numerically how much the approximation (|3.28p deviates from the precise error (j3.7p . and 
we will determine how fine the partial refinement needs to be. 



4. Numerical Method 

Based on the error estimator r], we need to decide what intervals at the qc-level shall be refined. 
To this end, we split r/ into individual contributions from each pc-level repatom: 



N-2 



where 



pc 

rjPj :-- 



ri= Yl C 

[gp" - jp'^ ip'impjpgp^]\jp'^ i'"i'^ j'^n^'^r'^ I'^'jiy"^)] , 

pc 



(4.1) 



(4.2) 



Next, we map the values rK to the qc-level intervals. Note that we already have rK = Q whenever 
J is a qc-repatom because we subtracted the qc-level interpolant from gP^ in (j3.28p . We accumulate 
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the remaining values ijj^ to the interval between the qc-repatoms j and j + 1 if the pc-repatom j 
lies between the qc-repatoms j and j + 1, that is, fij < j < /i^+i, to obtain 



qc 






(4.3) 



Moreover, we have to decide how to compute the partial refinement. Here we subdivide each qc- 
interval {lj,ij^i) into a fixed number A of subintervals of the same size, although other strategies 
are possible, like letting A depend on the respective interval size z^j = ij+i — ij, or dividing each 
interval into subintervals of different sizes. But here we stick to a fixed number A € N and an 
equidistant subdivision. We employ the following algorithm for partial refinement: 
j^-N + 1 

for j = -N + l,...,N -1 
to ^ max(l, I'j/A) 
0-1 ^0 

while ((72 < I'j) 

a\ <— min(o"i + uj, Uj) 

%^ ki-52 + ^J 

(T2 ^ CF2 + Z^J 

end 
end 
Note that this algorithm performs the necessary rounding if Vj is not divisible by A. If Uj < A, 
then the interval gets fully refined up to the atomistic level. 

We then use this partial refinement algorithm and the global and local error estimators above to 
construct the following algorithm for adaptive mesh coarsening: 

(1) Start with the model fully coarsened in the continuum region. 

(2) Solve the primal problem (I2.24p for y'^". 

(3) Determine the partial refinement pc according to the algorithm above. 

(4) Solve the dual problem ([3:25|) for gP". 

(5) Compute the global error estimator 77 from (|3.28p . 

(6) If r/ < Tgi, then stop. 

(7) Compute the local error estimators r/|'^ from (j4.2p and (|4.3p . 

(8) Refine all intervals {j,j + 1) with 



1 

- n 



rf:^ > max rfj^ (4.4) 



into two subintervals. 
(8) Go to (2). 

Here Tgi > denotes the global error tolerance. For the adaption criterion (14. 4p . numerical 
evidence has shown that Tfac = 10 is a reasonable choice. 

A possible improvement of step (8) would be to use the absolute value of the error estimator 
r/^'^ to decide into how many subintervals each interval {j,j + 1) shall be refined, instead of just 
refining into two subintervals. However, the magnitudes of the values r/^'^ differ considerably only 
in the first iterations. So this would only eliminate some of the early iterations which are still cheap 
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iteration 


#dof 


mill Uj 


max Uj 


ri 




|g(e-"'?-)| 


1 


12 


2048 


2048 


3.143618e-03 


3.143618e-03 


6.777614e-02 


2 


14 


1024 


1024 


5.208032e-03 


5.443530e-03 


6.463252e-02 


3 


16 


512 


1024 


8.771892e-03 


9.133002e-03 


5.946329e-02 


4 


18 


256 


1024 


1.293987e-02 


1.343519e-02 


5.074706e-02 


5 


20 


128 


1024 


1.520764e-02 


1.565599e-02 


3.787477e-02 


6 


22 


64 


1024 


1.267077e-02 


1.279361e-02 


2.271288e-02 


7 


24 


32 


1024 


6.760509e-03 


6.767672e-03 


1.004707e-02 


8 


26 


16 


1024 


2.395699e-03 


2.395699e-03 


3.286644e-03 


9 


28 


8 


1024 


6.933383e-04 


6.933394e-04 


9.216477e-04 


10 


32 


4 


1024 


2.061976e-04 


2.061988e-04 


2.638938e-04 


11 


40 


2 


1024 


5.841551e-05 


5.841804e-05 


6.391755e-05 


12 


54 


1 


1024 


7.567732e-06 


7.570401e-06 


9.376934e-06 



Table 1: Advance of the algorithm until the error toleranee Tgi — 10 ^ is achieved. 



due to the small number of unknowns involved there, whereas the computationally expensive later 
iterations will not be affected. Hence we expect the overall speedup to be small. 



5. Numerical Results 



Now we present and discuss our numerical results. We choose boundary conditions 
-M, yf| = -M + l, yl'i = M-l, y\\ = M. 






(5.1) 



The elastic moduli are given by fco = 0.1, k\ = 2, and ki = 1, and the lattice constant is oq = 1. 

We consider a chain of 4106 atoms, that is M = 2053. There is an atomistic region of 4 atoms 
from —1 to 2 around the dislocation at the center of the chain. The remaining part is modeled 
as continuum. We set atoms —3, —2, 3, 4 to be continuum repatoms so that ^'"^^(/"''y'?'^) and 
gaf^jaq^qc^ Can be evaluated without interpolation, and we set and atoms — M+1, — M+2, M—1, M 
to be continuum repatoms so that the boundary conditions can be set. Initially, the mesh in the 
continuum region is maximally coarsened, so that there are no more repatoms in addition to the 
ones already mentioned. This gives A^ = 6 and 12 qc-level degrees of freedom. There are two large 
elements of size v-^ = 1^4 = 2048, one on each side of the dislocation at the center of the chain, 
whereas the remaining elements necessarily have sizes I'j = 1 for j = — 5, — 3, — 2, — 1, 0, 1, 2,3, 5. 
The mesh is shown in the upper graph of Figure El 

The quantity of interest here is the size yi — yg of the dislocation at the center of the chain, 
that is, the distance between the two atoms and 1 to the left and right of the dislocation. The 
corresponding vector q E Vq reads as 



q= [0,..., 0,-1, 1,0,..., 0]^ 



(5.2) 



Table [T] shows how the algorithm proceeds for an error tolerance of Tgi = 10~^. One can easily 
see how the elements get refined and the error drops. Also shown is the precise error \Q{e°'^~'i^)\ 
which is available for this small model problem. Note that the column mmi'j refers only to the 
elements in the continuum region which actually can be coarsened, excluding the padding around 
the atomistic region and the boundary layer. 
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Figure 5: Mesh size Vj (green bar graph) and 77'*^ (black hnc graph) for error tolerances Tgi = 10 ^ (top), 
Tgi = 10^'"^ (middle) and Tgi = 10"^ (bottom). 



The bar graph m Figure [5] shows the size of the individual elements in the meshes before the first 
iteration and after the error tolerances Tgi = 10^^ and Tgi = 10^^ are achieved. One can clearly see 
how the elements get successively refined towards the center. 

The black lines in Figure [5] depict the decomposed error estimator 77^'^. As a result of the Galerkin 
orthogonality, it vanishes whereever the mesh is fully refined up to the atomistic level. We can as 
well read off the graphs how the algorithm tends to distribute the error uniformly over the whole 
mesh as the error tolerance is decreased. For Tgi = 10~^, we already achieved a close to flat region 
from elements —20 ... — 13 and 13 . . . 20. 
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g(gac-«c) 


A 


V 


EvT 


rj/ g(e— ^-) 


EvT/iQi^^'^m 


mesh 1 


6.777614e-02 


2 


3.143618e-03 


3.143618e-03 


0.046382 


0.046382 






4 


8.351650e-03 


8.351650e-03 


0.123224 


0.123224 






8 


1.708501e-02 


1.708501e-02 


0.252080 


0.252080 






oo 


6.777614e-02 


6.777614e-02 


1.000000 


1.000000 


mesh 2 


9.216477e-04 


2 


6.933383e-04 


6.933394e-04 


0.752281 


0.752282 






4 


8.749195e-04 


8.749264e-04 


0.949299 


0.949307 






8 


9.208978e-04 


9.209073e-04 


0.999186 


0.999197 






oo 


9.216477e-04 


9.216582e-04 


1.000000 


1.000011 


mesh 3 


9.376934e-06 


2 


7.567732e-06 


7.570401e-06 


0.807058 


0.807343 






4 


9.070422e-06 


9.085687e-06 


0.967312 


0.968940 






8 


9.320553e-06 


9.341074e-06 


0.993987 


0.996176 






oo 


9.376934e-06 


9.399414e-06 


1.000000 


1.002397 



Table 2: Efficiency of the error estimator, rj. 



Table [2] shows the efficiency of the error estimator for different meshes and different values of A. 
Meshes 1, 2, and 3 refer to the meshes displayed in Figure [5l For comparison, we also include the 
value A = oo, which indicates that the pc-level mesh for the dual solution is fully refined to the 
atomistic level. One can read off from this table that the error indicator gets closer to the precise 
error if the mesh gets finer. For the coarse mesh 1, the actual error is considerably underestimated. 
However, this has only little impact on the final mesh since all elements have to be refined anyhow. 
For the finer meshes 2 and 3, r/ gets closer to the precise error |(5(e"^~''^)|. 

Also, we can see from Table [2] how the choice of A affects the error estimator. As expected, ij 
gets closer to the precise error the larger A gets. For computational efficiency, we are interested 
in keeping A small, though. We can read off from the table for the fine mesh 3 that already the 
smallest possible value A = 2 gives a good estimate of the error, with a deviation of less than 20%. 
For A = 4, we already get a very precise estimate. 

The absolute accuracy of the error estimator is important, but even more important is how well 
r] controls the mesh refinement, that is how efficient the resulting mesh is in terms of reaching a 
prescribed error tolerance with a minimal number of degrees of freedom. We now investigate the 
influence of the parameter A on the mesh quality. 

Table[3]shows the number of degrees of freedom and the corresponding error during the automatic 
mesh adaption for different values of A. The choices A = 4, A = 8, and A = oo lead to different 
values of r/, but all result in the same mesh. For this reason, these values share a common column 
for T^dof and the precise error in the table. Figure [6] visualizes the relationship between the number 
of degrees of freedom and the error. We can see that A = 4, 8, oo leads to slightly better mesh that 
A = 2. However the difference is quite small so that the faster computation with A = 2 is very well 
acceptable. 
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Table 3: Mesh efBciency. 
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Figure 6: Mesh efficiency. 
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